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ABSTRACT 

We present a source and lens reconstruction for the optical Einstein ring gravitational 
lens system RXS J1131-1231. We resolve detail in the source, which is the host galaxy 
of a z = 0.658 quasar, down to a resolution of 0.045 arc seconds (this is the size 
of the smallest conclusively resolved structures, rather than the pixel scale), using 
a Bayesian technique with a realistic model for the prior information. The source 
reconstruction reveals a substantial amount of complex structure in the host galaxy, 
which is ^ 8 kpc in extent and contains several bright compact substructures, with the 
quasar source residing in one of these bright substructures. Additionally, we recover 
the mass distribution of the lensing galaxy, assuming a simply-parameterised model, 
using information from both the quasar images and the extended images. This allows 
a direct comparison of the amount of information about the lens that is provided by 
the quasar images in comparison to the extended images. In this system, we find that 
the extended images provide significantly more information about the lens than the 
quasar images alone, especially if we do not include prior constraints on the central 
position of the lens. 
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1 INTRODUCTION 

Gravitational lensing can be used as a powerful astrophysi- 
cal tool for probing the density profiles of galaxies, and is one 
of the few ways in which dark matter can be detected (e.g. 
iKoopmanslbOOSl ). In addition, it often magnifies source ob- 
jects by one to two orders of magnitude. This allows us to use 
the intervening gravitational lens as a kind of natural tele- 
scope, magnifying the source so that we can observe more de- 
tail than we would have been able to without the lens. This 
extra magnification provided by lensing has been very ben- 
eficial to studies of star formation and galaxy morphology 
at high redshift. Regions of the galaxy size and luminosity 
distribution that are inaccessible in unlensed observations 



are made (more) v i sible by lensing ( e.g. 'Pottini ct alj 2000l: 
Wavth et all l2005l: [Brewer fc Lewis! l2006b: Marsh all etHI 



20071 : iDve et all 120081^ The properties of the lens galax- 
ies (typically elliptical g alaxies) can also be inferred from 
their lensing effect (e.g. iKoopmans et al.l 120061 : iTreu et al.l 
120081 ). Of course, gravitational lensing distorts the image of 
the source, as well as magnifying it. Thus, techniques have 
been developed that aim to infer the mass profile of the lens 
galaxy and the surface brightness profile of the source, given 
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observ ed images (e.g. IWarren fc Dvell2003l : iBrewer fc Lewis! 
l2006a! ). 

The aim of this paper is to carry out this process with 
the recently discovered gravitatio nally lensed quasa r /host 
galaxy system RXS J1131-1231 l|Sluse et al.l |2003| ). This 
system consists of a quadruply imaged quasar at redshift 
z = 0.658 lensed by a galaxy ai z — 0.295. At the time of 
its discovery, it was the closest known lensed quasar, with 
some evidence for an extended optical Einstein ring - the 
image of the quasar host galaxy. Initial simple modelling 
suggested that the quasar source was magnified by a factor 
of ~ 50. Thus, subsequent ob servations with the AC S aboard 
the Hubble Space Telescope (|Claeskens et al.ll2006! . hereafter 
C06) allow the recovery of the morphology of the quasar's 
host galaxy down to a resolution of about 0.01 arc seconds 
- at least in principle, for the parts of the source nearest the 
caustic. Indeed, C06 presented a wide array of results based 
on HST observations (at 555nm and 814nm with ACS, and 
1600nm with NICMOS), including a detailed reconstruction 
of the extended source. 

The source reconstruction method used by C06 is based 
on lensing the image back to a pixellated grid in the source 
plane, setting the source surface brightnesses to equal the 
image surface brightness, and using a decision rule (in this 
case, the median) to decide on the value of a source pixel 
whenever two or more image pixels disagree about the value 
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of the same source pixel. If the point spread function (PSF) 
is small or the image has been deconvolved (in COG, the 
deconvolution step was neglected for the purposes of the ex- 
tended source reconstruction) and the lens model is correct, 
this method can expect to be quite accurate. However, in 
principle, the uncertainty in the lens model parameters and 
the deconvolution step should always be taken into account. 
In this paper, we focus our attention on the 555nm ACS im- 
age (the drizzled images, as reduced by C06, were provided 
to us), and the process by which we reconstruct the original, 
unlensed source from it. Any differences between our recon- 
struction and the C06 one can be attributed to the following 
advantages of our approach: PSF deconvolution, averaging 
over the lens parameter uncertainties, simultaneous fitting 
of all parameters, and the prior information that Bayesian 
methods are capable of taking into account: in the case of 
our model, that is the knowledge that the majority of pixels 
in an a strophysical sources should be dark iJBrewer fc LewisI 
l2006al ). The 555nm image is also of particular interest be- 
cause its rest wavelength (335nm) probes regions of recent 
star formation in a galaxy with an AGN. 

In the ca, s e of the Einstein Ring 0047-2808 
IjErewer fc LewisI l2006a ). our method was able to re- 
solve structure down to scales of ~ 0.01 arcsec, a factor of 
five smaller than that obtainable in an unlensed observation 
with the Hubble Spa ce Telescope and ab out double the 
resolution obtained bv lDve fc WarrenI ((2005) using adaptive 
pixellation and least squares applied to exactly the same 
data. This was possible because we used a prior distribution 
over possible sources that is more realistic as a model of 
our knowledge of an unknown astrophysical source, that is, 
we took into account the fact that it should be a positive 
structure against a dark background, a fact many methods 
(such as least squares an d some popular reg ularisation 
formulas) implicitly ignore (JBrewer fc LewisI l2006 al. These 
differences between methods are likely to be most significant 
when data are sparse or noisy, and all methods tend to 
agree as the data quality increases and we approach the 
regime where the observed image uniquely determines the 
correct reconstruction. 



2 BACKGROUND TO THE METHOD 

The conceptual basis of the Bayes i an reco nstruction method 
was presented in lBrewer fc LewisI (|200 6al. The idea is to fit 
a complex model to some data, but rather than simply op- 
timising the parameters of the model to achieve the best fit, 
we try to explore the whole volume of the parameter space 
that remains viable after taking into account the data. The 
effect of data is usually to reduce the volume of the plau- 
sible regions of the parameter space considerabljij. The ex- 
ploration of the parameter space can be achieved by using 
Markov Chain Monte Carlo (MCMC) algorithms, which are 
designed to produce a set of models sampled from the pos- 
terior distribution. In the case of modelling the background 
source and lens mass distribution of a gravitational lensing 
system, this allows us to obtain a sample of model sources 



^ For non-uniform probability distributions, "volume" is effec- 
tively the exponential of the information theory entropy of the 
distribution. 



and lenses that, when lensed and blurred by a PSF, match 
the observational data. The diversity of the models gives 
the uncertai nties in any q uantity of interest. The reader is 
referred to iGregorvl l|2005i ) for an introduction to Bayesian 
Inference and MCMC. 



3 METHOD AND ASSUMPTIONS 

The first step of a Bayesian analysis is to assign a likelihood 
function, or the probability density we would assign to the 
observed data if we knew the values of all of the parameters. 
To assign this, we need a noise frame, a measure of how un- 
certain we are about the noise level in each pixel. This is 
typically done by assuming that the observational error at 
pixel i is from a normal distribution with mean and known 
standard deviation ai. We extended this to include two "ex- 
tra noise parameters" a a and at, such th at the standard 
deviation for the error in the ith pixel is -v/cr^ + o"a + ""t/ij 
where fi is the predicted fiux in the j'th pixel, aa and at 
then become extra model parameters to be estimated from 
the data. The inclusion of aa and at implies that the extra 
noise level sigma varies with the predicted brightness of the 
pixel, with a square root dependence expected from Poisson 
photon counting. 

We chose the {ai} values to be zero for most of the 
image, but infinite for the brightest regions of the quasar 
images, effectively masking out those parts of the image; 
this mask can be seen in Figure [l] 

A model PSF was obtained using the Tiny Tim software 

(|Kristlll995 l'). As noted by C06, the TinyTim simulations did 
not successfully perform the geometric distortion correction, 
and the output PSF had slightly non-orthogonal diffraction 
spikes, whereas the spikes in the image are perpendicular. To 
correct this, the PSF was "straightened out" by evaluating 
it with respect to non- orthogonal axes; the resulting PSF is 
shown in Figure[2] Whilst this process is imperfect, the extra 
noise sigma protects against serious consequences resulting 
from slight inaccuracies in this process. While our choice 
of {ai} was designed to block out the brightest parts of the 
quasar images, since they are so bright, light from the quasar 
images still extends out past the masked regions and over- 
laps with interesting Einstein Ring structures. Thus, when 
modelling the image, we still required a flux component due 
to the quasar images. 

The four quasar images were modelled as being pro- 
portional to the corrected TinyTim profiles with unknown 
fiuxes and central positions. The surface brightness profile 
of the lens galaxy was modelled as the sum of two ellipti- 
cal Gaussian-like profiles (Sersic profiles, one for the core 
and one for the extended emission) proportional to e"*-*' , 
where R — ^^Qx''^ +y''^ /Q with unknown peak surface 
brightness, ellipticity Q, length scale L, and angle of ori- 
entation. The central position was also considered initially 
unknown, but for MCMC purposes the starting point was 
to have both profiles centred near the observed centre of the 
lens galaxy core. The slopes a were restricted to the range 
[0, 10] and assigned a uniform prior distribution, along with 
all of the other free parameters. Although elliptical galaxies 
are well modelled by a Sersic profile with a — 1/4, we are 
modelling this galaxy by two such profiles. This was done 
because the wings of the lens galaxy light profile (where it 
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overlaps the Einstein ring) are of great significance for our 
source reconstruction, and we do not want tlie core of the 
lens galaxy to be relevant to the wings. 

Note that there are parts of the image where flux is 
present from three sources: the lensed Einstein Ring, the 
wings of the PSF from the quasar images, and the fore- 
ground lens galaxy. The fact these all overlap suggests 
that the optimal approach (in all senses apart from CPU 
time) involves simultaneously fitting all of these compo- 
nents. Throughout this paper, any modelling has included 
all of the lens galaxy profile parameters as free parameter^ 
as well as the source, lens model parameters and positions 
and brightnesses of the four TinyTim PSF profiles, to model 
the contribution from the quasar images that remains even 
after masking out their central regions. 



4 LENS MODEL PARAMETERISATION 

The particular lens model we assumed for this sys- 
tem was a p scudo-isothcrmal elliptical potential (PIEP) 
IjBlandford fc Kochanck 1987), primarily for computational 
speed but also because it is fairly general and realistic, at 
least for single galaxy lenses that are not too elliptical. This 
model has five parameters: strength 6, ellipticity q (actually 
the axis ratio: g = 1 implies a circularly symmetric lens), 
orientation angle 6, and two parameters {xc, He) for the cen- 
tral position (measured relative to the central pixel of the 
images in Figure [!}. Although any Bayesian modelling can 
only explore a particular slice through the full hypothesis 
space we might have in our minds, using a simplified an- 
alytical model is often sufficient to illuminate the general 
properties of the true lens mass distribution. Also, it is typi- 
cally the case that inferences about the source of an Einstein 
Ring are insensi tive to the specific parameterisation for the 
lens model (e.g. IWavth et al.ll2005^ ■ as long as the model is 
able to fit the observed image at all. 

All Einstein rings can be expected to reside in an envi- 
ronment where the external shear due to neighbouring galax- 
ies is nonzero (Kochanek, private communication), and thus, 
external shear was also included in the lens model. There are 
two parameters for the e xternal shear : 7, its magnitude, and 
its orientation angle Oj. ISugai et al.l (|2007l ) have observed 
the flux ratios of the quasar images (via integral field spec- 
troscopy) and found that most of these ratios are consistent 
with a model of this type (elliptical potential plus external 
shear). A similar model was used by C06 (they used a sin- 
gular isothermal ellipse-l-7), where they find that it is the 
simplest parameterised model that can reproduce the obser- 
vations. In principle, we could adopt ever less restrictive pa- 
rameterisations for the lens, to hunt for substructures in its 
density profile. However, such an approach is extraordinarily 
computationally expensive (unless simplifying assumptions 
about the source are also made) and is beyond the scope of 
this paper. In terms of all of these parameters, the deflection 
angle formula, relating a point {x, y) in the image plane to 
a corresponding point {xs,ys) in the source plane, is 

Xs = X ~ ax{x,y) 



^ Except for the first part of Section \E\ where computational 
restrictions required that we fix the lens parameters. 
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Figure 2. Simulated PSFs from TinyTim. Each pixel corresponds 
to an ACS image pixel and is 0.049 arcseconds across. On the 
left is the actual PSF. For the purposes of the modelling of the 
extended images, this PSF was truncated to a 5x5 pixel PSF 
using only the brightest central parts of the PSF. For the quasar 
images, the wings are most important, and this can be seen most 
easily in the right hand panel. Since the quasars are so bright, 
the wings of the PSF are not negligible. 



Vs = y- 



Ax,y) 



(1) 



where the deflection angles a are given by the gradient of 
the potential 



i^{x,y) = -7(a;^ 



2- 

-y-Y. 



+ b^/qxfTyfJq 



(2) 



and (a;-y,j/-y) are the ray coordinates in the rotated coordi- 
nate system whose origin is (0, 0) and is aligned with the 
external shear (i.e. rotated by an angle O-y), and {xg,ye) are 
the ray coordinates in the rotated and translated coordi- 
nate system centred at {xc, 2/c) and oriented at an angle 0. 
The physical interpretation of each of these parameters sug- 
gests a plausible prior range for their values. To represent 
this knowledge we used the following prior distributions (Ta- 
ble [!}. Since these are broad distributions, and the data are 
good, the influence of these choices is negligible; they are 
included for completeness. 

In summary, the observed image was modelled as the 
sum of the following components: 

(i) Four TinyTim PSF profiles (Figure Ol with initially un- 
known amplitude and central position, to model the quasar 
images. While the bright parts of the images are masked 
out, the wings of the PSF are still important, so these com- 
ponents are required. 

(ii) Two elliptical Sersic profiles with initially unknown cen- 
tral position, peak surface brightness, scale radius, slope 
(Sersic index) and angle of orientation. One of these models 
the lens galaxy's core and the other models the fainter outer 
regions. 

(iii) A source, which is lensed by a PIEP-I-7 lens with un- 
known parameters and blurred by the 5x5 pixel core of the 
TinyTim PSF (Figure [U. In Section [01 the source is mod- 
elled in a simple way as a sum of six elliptical Sersic profiles 
(in order to find a good initial value for the lens parame- 
ters), and in section [6]the source is modelled as a pixellated 
grid with a prior distribution favouring non-negative pixel 
values and a dark background. 

(iv) Noise, to which we assign a Gaussian prior prob- 
ability distribution with unknown standard deviation 
^/o■f + ai + a-ffi for each pixel, with aa and ab initially 
unknown and the ai specifled in advance to mask out the 
pixels with the most systematics. 
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Figure 1. The logarithm of the observed image (on the left) shows both the quasar images, along with their additional effects such 
as diffraction spikes, and the faint Einstein ring image of the host galaxy. On the right is the image (scaled linearly) with some parts 
blocked out - these blanked regions are those where the {cri} have been adjusted to block out the brightest emission from the quasar. 
The (Ti for some pixels has also been artificially boosted to reduce the effects of the outer parts of the diffraction spikes from the quasar 
images. We expect the inner parts of these spikes to be well modelled by the TinyTim profiles (Figure [2J. Note that the lens galaxy light 
profile extends over this entire image; where it appears that the image is blank, there is actually a positive fiux. 

Table 1. Prior probability densities for the lens model parameters, and also the extra noise parameters era and ct;,. 



Parameter 



Prior Distribution 



b Normal, mean f.8", SD 0.5, 6 > f 

q Normal, mean 0.9, SD 0.2, < g < 1 

Xc Normal, centred at the lens galaxy core, SD f .0" 

j/c Normal, centred at the lens galaxy core, SD f .0" 

6 Uniform, between and it 

Ca Improper Uniform, Ca > 

(T(, Improper Uniform, (Tj, > 

0.y Uniform, between and 2it 

J Exponential, mean 0.1 



One of the most difficult tasks in a source reconstruc- 
tion problem is to find good values of the lens parameters to 
serve as the starting point for an MCMC run. If the source 
is pixellated, then exploration of the lens parameter space 
is slow because we effectively have to marginalize over thou- 
sands of source dimensions - so if we start with incorrect 
lens model values, the burn-in approach to the more plau- 
sible values will be slow. A good starting point for the lens 
can usually be obtained by running a much simpler version 
of the whole inference problem, for instance with a simpler 
model for the source, or by using only the quasar image po- 
sitions and brightnesses as constraints. In the next section, 
we apply the latter approach to see how well the quasar im- 
ages alone can constrain the lens model parameters. Later 
(Section 15. 3p , the extended images are taken into account 
by using a simply parameterised source model, where the 
number of source dimensions to be marginalized is only 36. 
Finally (Section |6}, we use a pixellated model for the source 
in order to reconstruct its structure with minimal assump- 
tions. 



5 HOW THE QUASAR IMAGES CONSTRAIN 

THE LENS 

5.1 Theory 

The four quasar images can constrain the lens model be- 
cause we require that the four image positions lens back to 
the same point in the source plane. Actually, with an up- 
per limit of ~ 0.02 pc for the co ntinuum source size (e.g. 
IWavth. O'Dowd. fc Websteij |200a ) . this exact requirement 
is too strong - we can really only insist that they lens back 
to within ~ 3 microarcseconds of each other (using a con- 
cordance c osmology Qm = .27, f^A ~ 0.73, Ho ~ 71 km 
s"^Mpc~\ iNoha et all (|2004l )). The results of this subsec- 
tion are unchanged i f we use a smaller quasar size of lO'^^pc 
(jMorgan et al.l 120081 ) because the limiting factor is the ac- 
curacy of the astrometry, rather than the assumed quasar 
source size. The magnificatio ns of the images ca n also pro- 
vide some information (e.g. iLewis et al.l ( 20021 ) used the 
brightness of the third image of a quasar to argue for a 
model for the lens mass profile that creates a naked cusp in 
the source plane) although microlensing and absorption ef- 
fects can limit the usefulness of including the magnifications 
as constraints for more typical systems. 



Lens Inversion of RXS J 113 1-1231 5 



The quasar positions were measured by fitting Gaus- 
sians to the peaks of the images. For the purposes of cen- 
troiding, this is an adequate approximation: the alternative 
is to calculate very high resolution simulated PSFs. Differ- 
ent (unknown) noise levels were assumed for each of the four 
images; however, the uncertainties in position were found to 
be ^ 0.003 arc seconds for all of the four images. This cor- 
responds to less than 1/lOth of an image pixel. 

To infer our PIEP-I-7 parameters from the quasar po- 
sitions, we first implemented an MCMC algorithm to ex- 
plore the prior distribution for the lens parameters (Ta- 
ble [T|. Then we imposed a constraint on the probabil- 
ity distribution: that the expected value of the mean 
inter-pair distance, or spread, of the four points in the 
source plane upon back ray tracing, should be about 
10~^ arc seconds. This modifies the probability distribu- 
tion over the lens parameter space by a multiplicative factor 
exp( — fc X Spread(Lens Parameters)), where the value of 
k is chosen so that the expectation value takes the value 
we wish to impose, 10~^ arc seconds. Conventionally, one 
would estimate the lens parameters by finding those values 
that minimise the spread in the source plane. Our proba- 
bilistic approach softens this constraint and implies that we 
sample from the range of lens models that reduce the scatter 
to about 10~^ arc seconds or less. 

This approach assumes that we know the true exact 
image positions, although it can be extended to allow for 
uncertainty in the image positions, as follows. Denote the 
true unknown quasar image positions collectively by X, the 
estimated positions by x and the lens parameter values by L. 
By the rules of probability theory, the posterior probability 
distribution for L given x (here, we assume that the known 
data are the a;'s, rather than the entire image) can be written 
as: 



p{L\x) 



p{L,X\x)dX 



p{X\x)p{L\X,x)dX 



(3) 



(X / p(X|2:)p(i)exp(-fc X Spread(X, L))dX 



Since knowledge of X would make x irrelevant, this becomes 
p{L\x) = / p{X\x)p{L\X)dX 

(4) 

Hence, the true image positions {X} can be introduced as 
extra nuisance parameters to be estimated, and then we can 
sample the distribution under the integral sign in Equa- 
tion m The small ccntroiding uncertainties of ± ~ 0.003 
arcseconds were taken to specify the standard deviations of 
Gaussian probability distributions for p{X\x). 

The reader may wonder whether it would be more cor- 
rect to introduce the unknown source plane position of the 
quasar as the nuisance parameters, calculate the image po- 
sitions using the lens model, and use the x's and their error 
bars to define a likelihood function. Whilst there is nothing 
wrong with this approach, it does involve the computation- 
ally challenging task of inverting the lens equations (Equa- 
tion [1]) to find the image positions. Our source plane spread 
approach is much easier to implement computationally, but 



relies on the unconventional step of directly assigning a pos- 
terior probability distribution: p{X\x). 

5.2 Results 

To implement the inference described in the previous sec- 
tion, a Metropolis-Hastings sampler was written to tar- 
get the posterior distribution of Equation |4l Unfortunately, 
this simple implementation had serious drawbacks. The 
joint posterior distribution for the lens parameters and true 
quasar image positions consists of long, thin, curved tun- 
nels of high probability, and most standard sampling tech- 
niques have very poor mixing properties when sampling from 
such highly correlated distributions. To overcome this chal- 
lenge, we used a different sampling t echnique, Linked Im- 
portance Sampling (LIS) (JNeaJ [2005|) . LIS produces inde- 
pendent weighted samples from the target probability dis- 
tribution. It only requires that we can independently sample 
from a simple distribution (e.g. the prior) and can define 
valid MCMC transitions with respect to distributions that 
are in some sense intermediate between the prior and poste- 
rior. For example, the common 'annealing' approach of rais- 
ing the likelihood to some power (< 1) can be used within 
the LIS framework. The fact that each LIS run produces 
an independent sample from the target distribution makes 
it an attractive option for sampling from highly correlated 
distributions. The only possible pitfall is that the weights 
can vary significantly, such that the sample is completely 
dominated by only one or two points with large weight. 

We repeated this analysis twice: first, with the conser- 
vative priors on the central position of the lens - the priors 
in Table [1] For comparison, we ran the algorithm with much 
more informative priors on the central position {xc, yc) of the 
lens, an uncertainty of 1 pixel or ±0.049 arcseconds. The re- 
sults are shown as the blue and red curves in Figure [3] Far 
from uniquely determining the lens model, the quasar images 
have only managed to give moderately strong constraints on 
the overall strength of the lens (6) and the angle of orienta- 
tion of the external shear {6~f). For all other parameters, the 
marginal distributions are very wide, in some cases nearly as 
wide as the priors, so the quasars have provided only a small 
amount of information about them. In the seven dimensional 
space of the lens parameters, the quasar data yields a pos- 
terior distribution that contains long, narrow tunnels: while 
the volume of possible lens models is significantly decreased, 
the degeneracies inherent in lensing prevent precise inference 
about any single parameter. 

With the stronger prior information about the central 
position (red curves), the marginal probability distributions 
tighten substantially, but not well enough to give a reliable 
lens parameter estimate for which we could trust a source 
reconstruction. In the next section, we use a simplistic model 
for the extended source in an attempt to get a good starting 
estimate for the lens model parameters. This estimate will 
then be used in the final run (Section [6]) with a pixellated 
source plane. 



5.3 Simply- Parameterised Modelling of the 
Extended Source 

Rather than using pixels from the outset, we first modelled 
the source as the sum of six elliptical Gaussian-like (Ser- 
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Figure 3. Comparison of the inferred lens parameters (estimated marginal probability densities, unnormalized) based on three data 
sets: (1) The quasar image positions and uncertainties, with a weak prior distribution for the central position of the lens (blue curves). 

(2) The quasar image positions and uncertainties, with a strong prior (it 1 pixel) on the central position of the lens (red curves), and 

(3) The entire image, and a pixellated source model, run at a temperature T=10 (black curves). Note that the parameter spaces for the 
two angles 9 and 9-y are periodic with period it. 



sic) profiles, with unknown brightness, orientation, central 
position, slope and ellipticity for each. Similarly, the lensing 
galaxy light profile was modelled as two elliptical Sersic pro- 
files in the image plane. This simplified model reduces the 
dimensionality of the problem from thousands to 36, and 
makes it much more likely that a simple search algorithm 
can find something close to the optimal values for the lens 
model parameters. We used a simple Metropolis algorithm 
to derive estimates and uncertainties on the lens parameters. 
This estimate was used as a starting point for the chains in 
Section [6l where we use a pixellated source. 

A typical simple source model from the sample is shown 
in Figure |4l Within this parameterisation, the uncertainties 
about the source are very small and Figure |3] can essen- 
tially be interpreted as the unique source reconstruction. 
The scales on the axes are the same for the source plane 



and the image plane, so an idea of the magnification can 
be obtained visually. Given this model, the data favour a 
highly complex source (since the blobs do not overlap to 
the extent that they become indistinguishable), lensed by 
a slightly elliptical lens whose centre is located close to the 
centre of the observed lensing galaxy. The individual Sersic 
profiles in the source shown in Figure U have been colour 
coded and labelled, making identification of the correspond- 
ing images easier. These labels will be used throughout the 
paper to refer to specific substructures in the source, and 
their corresponding images. 

Comparing the image in Figure |4] to the observed one 
in the right hand panel of Figure [l] we see that all of the 
basic structures observed have been reproduced by this sim- 
ple model. However, the simply-parameterised model cannot 
reproduce the exact shapes of the observed images. For ex- 
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ample, the part of the source labehed A should contain sub- 
structure, because we can see that the simple model has pre- 
dicted a continuous image Al, yet the actual data contains 
blank gaps in some parts of that image. Another example is 
that image B2 has irregular brightness variations along its 
length, something the simply parameterised source model 
cannot reproduce. There are also some very faint additional 
structures that have been missed, such as a third faint inner 
ring below A2, that could be a continuation of the image 
Dl. These differences can also be seen in the residuals in 
Figure (5] 

The bright ring (images CI and C2) passes through 
the quasar image positions, so the quasar source is located 
inside source component C, but within the diamond caustic 
since the quasar has been imaged 4 times. However, due to 
uncertainty in the lens parameters, the quasar source cannot 
be located more accurately than this in the source plane. 
Source component C is moderately elliptical and is about 
0.2 arcseconds in length, corresponding to a physical length 
scale of ~ 1400 pc. The estimated magnifications for each 
component of the source are as follows: A (12.3), B (12.2), 
C (20.0), D (3.9), E (7.5), F (4.8), Quasar (45.0). 



only change in extremely tiny steps: for example, typically 
by 10~^, as far as the data and the current source will al- 
low, so the lens model would explore the marginal distri- 
bution for lens parameters very slowly. Unfortunately, LIS 
was found to be unfeasible (the weights varied by orders 
of magnitude) due to the massive number of source pa- 
rameters used. Hence, there was no realistic option other 
than to fix the lens model at a reasonable value and run 
the MCMC for the source variables only. We fixed the lens 
at the best values obtained from the simply-parameterised 
exte nded source model, and then r an a slow "annealing" 
( Kir kpatrick. Gelatt. fc Vecchil [l983J) simulation to deter- 
mine a good lens model, which was then fixed, and the source 
space was explored at an annealing temperature of 1, to sam- 
ple from the posterior distribution for the source pixel values 
given the lens. A large number of simulations with slightly 
different lens parameter values was done to verify that the 
source reconstruction is insensitive to slight uncertainties in 
the lens parameters. The results are presented in the next 
section. 



6 PIXELLATED SOURCE MODEL 

To obtain high resolution imaging of the source, we divided 
the source plane into a grid of 200x200 pixels, each having 
a width of 0.015 arc seconds. The MCM C algorithm used 
is almo st the same as that described in iBrewer fc Lewis! 
l|2006bl ). and is described in greater detail there. Starting 
from a blank source, proposal changes are made that either 
add a bright "atom" of light to the source plane. Each atom 
has a position in the source plane indicating which pixel it 
is in, as well as a positive brightness, with exponential prior 
distribution. Small proposal transitions can be made which 
slightly adjust the parameters of an atom, in accordance 
with our chosen prior for these parameters. If the proposed 
source is a better fit, it is accepted, if it is a worse fit it is 
accepted with an acceptance probability equal to the ratio 
of the proposed likelihood to the curre nt likelihood; t his is 
just the standard Metropolis algorithm ([GrcEorr'2005) with 
a Massive Inference prior distribution (Skilling 1998) for the 
source. 

The prior expected value for an atom's flux can either 
be specified in advance, or incorporated as an additional 
hyperparameter to be estimated. The latter approach is at- 
tractive and has been used in the context of Gaussian priors 
for the source pixel values (|Suvu et al.ll2006h . We chose the 
former approach for increased computational efficiency, and 
found that different values for this hyperparameter did not 
significantly change the final reconstruction, provided that 
the value was not so low that the reconstructed sky was 
bright or so high that the model cannot detect the presence 
of structures that are obviously real. 

A straightforward implementation of any MCMC 
method is highly inefficient for the problem we have just 
posed. This is because we need to be able to explore the 
marginal posterior distribution (marginalising over possi- 
bly thousands of source parameters) for the lens parame- 
ters with reasonable efficiency. If we just alternated between 
source updates and lens updates, the lens parameters would 



7 RESULTS 



T.l Source 



A final estimate for the source is shown in Figure [SI which is 
obtained by averaging all sources encountered by the Markov 
chain. Although the diversity of the samples encountered is 
an accurate representation of the level of uncertainty about 
any pixel values, it is inconvenient to present a large sam- 
ple of images. Additionally, the spiky fluctuations caused by 
the Massive Inference prior remain in the sampled sources. 
Taking the mean of all of the sources provides a single es- 
timate of the source profile that is optimal in the sense of 
minimising the expected square error; it also creates a more 
visually appealing smooth reconstruction where all of the 
spiky fluctuations in the individual samples have been aver- 
aged out. It is also quite natural for uncertain areas to be 
blurred in images, and for people viewing images to interpret 
a smooth image as possibly being caused by an underlying 
complex image. 

The positions of the major bright substructures in the 
source (Figure[6)l are in agreement with those found by C06. 
However, our reconstruction presents a clearer view of the 
compact central source C, where the quasar resides. In the 
simply parameterised case, source C was found be elliptical. 
In the pixellated reconstruction, this part of the source has 
an elliptical component but with extra nearby structures 
that follow the caustic. The images of this extra structure lie 
within images CI and C2. This suggests that the elliptical 
component is sufficient to explain the position of the ring 
C2 but not its brightness; the algorithm tries to account for 
this by adding extra source flux in regions that have images 
within C2 but only within those parts of CI that have been 
masked out. It seems more plausible that image CI is more 
affected by dust absorption than C2, rather than the source 
coincidentally following the caustic. In principle, we could 
parameterise the unknown dust profile of the lens galaxy 
and simultaneo u sly es timate this from the data (as done by 
e.g. ISuvu et aU l|2008l )'). but this is beyond the scope of the 
present paper. 
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Extended Source 



Lensed, Blurred Image 





Figure 4. Reconstruction of the extended source with six eUiptical Gaussian-like Sersic profiles, shown on the same scale as the observed 
image. Most of the basic features of the observed image have been reproduced, but not their detailed shapes. The differences tell us where 
we should expect to see some additional substructure in the pixellated source model (see text). The positions of the caustics and critical 
curves are shown in white, although there should really be a small amount of uncertainty about their positions due to the uncertainty 
in the lens parameters. 



Predicted Image {including lens galaxy) 



Standardised Residuals 




-2-1 1 2 
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Figure 5. Model-predicted lensed, blurred image, using a typical model sampled from the posterior distribution, and the simply- 
parameterised source model. On the right are the standardised residuals. While the basic features of the observed image are reproduced, 
there are details that have not been well modelled by the simply-parameterised source. 



The predicted gap in the middle of source component 
A is also clearly present in our reconstruction, whereas it 
is much less clear, if present at all, in the reconstruction by 
COG. Source component E also contains some substructure of 
its own. In addition, the surface brightness contrast between 
these features and the diffuse background is much greater in 
our reconstruction, although this is simply because we are 
focusing on the 555nm image. At other wavelengths, the 
compact substructures are less pronounced. 

To reduce these systematic effects, we repeated the sim- 
ulations at an annealing temperature of 10, to allow freer 
exploration of the parameter space. This is an ad-hoc de- 
vice that is not as justified as explicit modelling of dust 
(|Suvu et al.ll2008r ). but is helpful nonetheless. At a temper- 



ature of 10, it also becomes computationally feasible to free 
the lens parameters. The results, shown in Figure |S1 still 
show clear images of the bright, compact substructures that 
are present in the source. These bright substructures account 
for all of the bright images that are visible in the data; the 
diffuse emission that the T=10 analysis misses is caused by a 
thicker, fainter ring (image E2) that is not the most striking 
feature of the image. 

The reconstruction produced by the high temperature 
run still contains source plane structure that follows the 
caustic. This suggests that increasing the temperature has 
not completely negated the effect of dust. However, the 
raised temperature simulation still results in model images 
that reproduce the positions and shapes of all of the details 



Lens Inversion of RXS J 113 1-1231 9 
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Figure 6. Posterior mean source, with fixed lens parameters. On the right is a greyscale version of the original reconstruction by C06, 
for comparison. Note that the C06 reconstruction used data in three filters, which partially accounts for its more diffuse appearance: the 
compact substructures are most notable in the 555nm image (the subject of this paper). 
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Figure 7. Model-predicted lensed, blurred image, using a typical model sampled from the posterior distribution, and the pixellatod 
source model. On the right are the standardised residuals. 



in the observed Einstein ring, while being more permissive 
about their fluxes. Thus, the source reconstruction in Fig- 
ure [S] is a good model for the positions and (excluding the 
central source C) shapes of the source galaxy substructures. 

This system is probably one of the most complex Ein- 
stein Rings known, with a spectacular number of distinct 
extended images. Hence, it is unsurprising that the source 
galaxy morphology is also very complex. The redshift of 
the source {z = 0.658) implies that this hght was emitted 
in the near UV (wavelength of 335 nm), suggesting that 
these structures are mapping the star forming regions of the 



source galaxy. With the assumed cosmology, the length scale 
in the source plane is ~ 6.959 kpc/arcsecond. The bulk of 
the galaxy is just over 1 arc second across, so the entire 
source is about 8 kpc across; the source is a medium-sized 
galaxy. Compared to C06, our reconstruction is devoid of a 
large amount of extended emission. This is probably caused 
by some Ught from the lens galaxy being attributed to the 
source in their analysis. The surface brightness ratio of the 
compact structures to the diffuse emission can be estimated 
visually by simply looking at the image; where the contrast 
is a lot stronger than it appears to be in the COG reconstruc- 
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Temperature = 10 Mean Source 
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Figure 8. Mean source encountered by a chain run at a temperature of 10. Tliis allows tlie cliain to explore the parameter space more 
freely, reducing the effect of unmodelled systematics present in the Temperature=l reconstruction (albeit by discarding information). 
Most of the substructures are still sharply resolved. On the right, the posterior probability that any pixel is nonzero is plotted. White 
corresponds to 0, while black corresponds to 1. While it is very difficult to be confident about any particular pixel (due to the finite 
information content of the image), groups of pixels with consistently high values of this probability represent a secure detection. 



tion. There is also evidence for a companion dwarf galaxy 
2.4 kpc away and roughly 700 pc in diameter, to the left of 
the main galaxy (source F). We have also found that the 
quasar resides within an elliptical {q ~ 0.5) region of young 
stars that is about 1400 pc in extent. 



7.2 Lens Parameters 

When running the MCMC simulation at a temperature of 
10, inferring the lens parameter uncertainties becomes com- 
putationally feasible. Thus, we can compare the lens con- 
straints derived from the extended images with those derived 
from the point-like images. The marginal posterior distribu- 
tion for the lens parameters, with a pixellated source, are 
shown as the black curves in Figure [3] This shows clearly 
that the extended images provide significantly stronger con- 
straints on the lens model parameters in this system, in 
contrast to the claim made by COG that the opposite is 
true, and in agreement with iKochanek. Keeton. fc McLeodI 
(|2001f ). Note that the distributions for h inferred from the 
quasar astrometry and the extended image reconstruction 
overlap only slightly. This is because the very small quoted 
astrometric uncertainties of 0.003" (less than 1/lOth of an 
image pixel) do not take into account known systematic 
effects such as the presence of background flux from the 
lens galaxy and the Einstein Ring, the fact that the peak 
of the PSF is not really a Gaussian, and the fact that the 
lens is not really a PIEP-I-7. This small disagreement means 
that the parameters inferred from the extended images are 
only marginally capable of giving quasar image positions 
correctly to within 0.003". They are, however, capable of 
reproducing the positions to within a relaxed tolerance of 
~ 1/3 of a pixel. The axis ratio q, describing the ellipticity 
of the lens potential, is found to be 0.935 ± 0.005. Clearly, 
the data rule out q = 1 and therefore Singular Isothermal 
Sphere -I- 7 models (this is still the case even when only the 
quasar data is used). If we had assigned a delta function 
of prior probability at g = 1, the data would downgrade 



its posterior probability significantly when compared to any 
realistic diffuse prior for q. 

The mass of the lens can be calculated from these re- 
sults. For the purposes of this calculation, the lens is ap- 
proximated as circular. For an isothermal sphere, the de- 
flection angles at the Einstein radius of the ring (now equal 
to 6) is simply b. For a point mass at the origin, the re- 
quired mass to produce the same deflection would be h^ . 
Since lensing obeys Gauss' law, 6^ must also be the amount 
of mass enclosed within the ring. An alternative approach 
would be to calculate the nondimensional density, which 
proportional to the 2-D Laplac ian of the lens potential 
(jSchneider. Ehlers. fc Falcd 119921 ). In scaled units, the to- 
tal mass of the lens is therefore h^ — 3.205 ± 0.01, where the 
uncertainty was found, as usual, by considering an ensem- 
ble of lens models from the MCMC output and calculating 
the mass for each. However, systematics introduced by the 
approximations in computing this value will probably over- 
whelm this quoted uncertainty. The mass unit that would 
give an Einstein Radius of 1 arc second needs to be com- 
puted to convert this figure into physical units. When this 
is done, the estimated mass of the lensing galaxy (within the 
ring; the total mass outside of this is very poorly contrained 
by lensing) is found to be (6.95 ±0.02) x lO" Mq. With the 
isothermal assumption, the velocity dispersion is ~ 350 km 
s"'^. From the lens galaxy Hght profile parameters, we find 
that the flux of the lens galaxy (within the Einstein Ring) is 
close to 50% of the flux of the brightest QSO image, which 
has a magnitude of 17.74 in the 555nm fllter. At a redshift 
of 2=0.295, the luminosity distance to the lens (with our as- 
sumed cosmology) is 1510 Mpc. Thus, the average (within 
the ring) mass to Hght ratio of the lens galaxy is found to 
be 8.8 Mo/Lo. 

The lens potential is only slightly elliptical and is cen- 
tred near the centre of the lens galaxy, if the centre is deflned 
by the brightest pixel of the galaxy core. All of these con- 
clusions are similar to those made of ER 0047-2808, and are 
probably typical features of Einstein rings and all other sys- 
tems with single galaxy lenses. This lends more support to 
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the often used assumption that the centre of a lens galaxies 
light profile is also the point at which the lens model should 
be centred. It remains unclear how the galaxy light profile 
information could be taken into account in a more complex 
kiloparametric model of the lens; that is a topic for further 
research. The lens light profile parameters are presented in 
Tabled 

A preliminary investigation of the time delays predicted 
by our lens model suggests that it does not exa c tly rep roduce 
the time delays measured by [Morgan et al.l l|2006l'). Since 
this system exhibits sig nificant microlensing (|Chartas et al.l 
120081 : ISluse et al.l 120081 ) , the time delay measurements are 
uncertain, but it is possible that the PIEP+7 lens model 
will be ruled out by further observations. This would not 
be catastrophic for the present study, for several reasons. 
Firstly, source reconstructions tend to b e insensitive to sligh t 
misspecification of the lens model (e.g. IWavth et al.l 12005!). 
Secondly, all parameterisations are false. We already know 
from prior information that the lens is not really a PIEP+7 
model. All modelling can only consider a single "slice" 
through a full hypothesis space, and the conclusions reached 
on that slice may or may not be representative of the full 
space. They often are, but there are never any guarantees. 



constrain the lens model. Simultaneous multi - wavelength re- 
const ructions are now becoming routine (e.g. [Marshall et al.l 
[20071). However, all of the structures in these images are in 
the same locations, and so a multi-wavelength reconstruc- 
tion would not produce significantly different conclusions to 
those reached here. COG note that in the near infrared image, 
the compact bright images are less pronounced compared to 
the diffuse background, which is what would be expected if 
the substructures are regions dense in hot young stars. 

This study has relied on a immber of common as- 
sumptions that future research will seek to relax. Extending 
lens reconstruction techniques to incorporate kiloparametric 
models of the source and the lens simultaneously is an am- 
bitious ta sk, but some step s are already being taken in that 
direction IjSuvu et al.ll2008r ). Flexible lens modelling plus in- 
formation from time delay measurements and other sources 
would be extremely valuable for studies of galaxy dark mat- 
ter haloes. Also, explicit modelling of dust absorption by the 
lens galaxy is proving to be an important ingredient in the 
inversion of Einstein Rings and would be an essential part 
of future work on this system. 



8 CONCLUSIONS AND FURTHER WORK 

In this paper, we have presented a detailed gravitational 
lens reconstruction of the optical extended source in the 
Einstein Ring RXS J1131-1231. The source is a medium 
sized galaxy {^^ 8 kpc in visible extent) with several compact 
bright emission regions. The substructures we found are in 
general agreement with those found by COG in terms of their 
position, but we have shown that they are brighter and more 
compact than was previously determined. In addition, our 
reconstruction provides a clearer view of the substructures, 
including near the central regions of the source. The quasar 
resides in a bright emission region with an extent of about ~ 
0.15 arcseconds or 1 kpc. It should be noted that the wave- 
length of the observations in the rest frame is 335 nm, so 
this reconstruction traces regions of recent star formation in 
the source galaxy. 

We have also directly compared point images vs ex- 
tended images with regard to how well each is able to con- 
strain the lens model. We found that there is a significant 
gain to be made in taking into account all of the informa- 
tion from the extended images. It has been sug gested that 
this is not true in general (Fcr reras. Saha. fc B uries 2007)), 
although it really depends on the resolution and number 
of extended images, which in this case is high. Certainly, 
in using both, there is nothing to lose but CPU cycles. 
This system has the potential to become one of the most 
well-constrained gravitational lenses, with multiple images 
of the extended ring, quasar image positions and flux ratios 
in multiple bands, and time delay measuremen ts available 
([Morgan et al.ll2006l : JKeeton fc Moustakad 120051 ). Hence, it 
should be possible to carry out a detailed kiloarametric 
study of its mass proflle to shed some light on the dark 
matter halo of the lens galaxy. 

This paper was based on a single image of this system, 
the 555nm ACS image. Other HST images at different wave- 
lengths (814nm, l.G^m) are available (COG) and can further 
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